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A nonparametric Bayesian approach is developed to determine quantum potentials from empirical 
data for quantum systems at finite temperature. The approach combines the likelihood model 
of quantum mechanics with a priori information over potentials implemented in form of stochastic 
processes. Its specific advantages are the possibilities to deal with heterogeneous data and to express 
a priori information explicitly, i.e., directly in terms of the potential of interest. A numerical solution 
in maximum a posteriori approximation was feasible for one-dimensional problems. Using correct 
a priori information turned out to be essential. 
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The last years have seen a rapidly growing interest 
in learning from empirical data. Increasing computa- 
tional resources enabled successful applications of empiri- 
cal learning algorithms in many different areas including, 
for example, time series prediction, image reconstruction, 
speech recognition, and many more regression, classifica- 
tion, and density estimation problems. Empirical learn- 
ing, i.e., the problem of finding underlying general laws 
from observations, represents a typical inverse problem 
and is usually ill-posed in the sense of Hadamard . 
It is well known that a successful solution of such prob- 
lems requires additional a priori information. In the set- 
ting of empirical learning it is a priori information which 
controls the generalization ability of a learning system by 
providing the link between available empirical "training" 
data and unknown outcome in future "test" situations. 

The empirical learning problem we study in this Pa- 
per is the reconstruction of potentials from measuring 
quantum systems at finite temperature, i.e., the prob- 
lem of inverse quantum statistics. Two classical research 
fields dealing with the determination of potentials are 
inverse scattering theory Q and inverse spectral theory 
1^,^ . They characterize the kind of data which are nec- 
essary, in addition to a given spectrum, to identify a po- 
tential uniquely. For example, such data can be a sec- 
ond complete spectrum for different boundary conditions, 
knowledge of the potential on a half interval, or the phase 
shifts as a function of energy. However, neither a com- 
plete spectrum nor values of potentials or phase shifts 
for all energies can be determined empirically by a finite 
number of measurements. Hence, any practical algorithm 
for reconstructing potentials from data must rely on ad- 
ditional a priori assumptions, if not explicitly then im- 
plicitly. Furthermore, besides energy, other observables 
like particle coordinates or momenta may have been mea- 
sured for a quantum system. Therefore, the approach we 
study in this Paper is designed to deal with arbitrary, 
also non-spectral, data and to treat situation specific a 
priori information in a flexible and explicit manner. 

Many disciplines contributed empirical learning al- 
gorithms some of the most widely spread being de- 



cision trees, neural networks, projection pursuit tech- 
niques, various spline methods, regularization ap- 
proaches, graphical models, support vector machines, 
and, becoming especially popular recently, nonparamet- 
ric Bayesian methods . Motivated by the clear and 
general framework it provides, the approach we will rely 
on is that of Bayesian statistics |p^-p^ which can easily 
be adapted to inverse quantum statistics. Computation- 
ally, however, its application to quantum systems turns 
out to be more demanding than, for example, typical ap- 
plications to regression problems. 

A Bayesian approach is based on two probability den- 
sities: 1. a likelihood model p{x\0,v), quantifying the 
probability of outcome x when measuring observable O 
given a (not directly observable) potential v and 2. a 
prior density po{v) — p{v\Do) defined over a space V 
of possible potentials assuming a priori information Dq. 
Further, let Dt = {xt,Ot) = {{xi,Oi)\l < i < n} denote 
available training data consisting of outcome-observable 
pairs and D = Dt U Dq the union of training data and 
a priori information. To make predictions we aim at cal- 
culating the predictive density for given data D 

p{x\0,D) ^ Jdvp{x\0,v)p{v\D). (1) 

The posterior density p{v\D) appearing in that for- 
mula is related to prior density and likelihood through 
Bayes' theorem p{v\D) = p{xt\Ot,v) po{v)/po{xt\Ot), 
where the likelihood factorizes for independent data 
p{xt\Ot,v) = YiiPi^il^ij''^) denominator is v- 

independent. The w-integral stands for an integral over 
parameters if we choose a parameterized space V, or for 
a functional integration if V is an infinite dimensional 
function space. Because we will in most cases not be 
able to solve that integral exactly we treat it in max- 
imum a posteriori approximation, i.e., in saddle point 
approximation if the likelihood varies slowly at the sta- 
tionary point. Then, assuming p{x\0,D) « p{x\0,v*) 
with V* — argmax^gi;p(ti|£') integration is replaced by 
maximization. 

According to the axioms of quantum mechanics observ- 
ables O are represented by hermitian operators and the 
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probability of finding outcome x measuring observable O 
is given by 

pix\0,v)=TT(^Po{x)piv)y (2) 

Here p denotes the density operator of the quantum me- 
chanical system and Poi^) = X^iy l^(i^) is the 
projector into the space spanned by the eigenfunctions 
> of O having eigenvalue x and which we assume 
to be orthonormalized. If a system is not in an eigen- 
state of the observable, a quantum mechanical measure- 
ment changes the state of the system. Hence, to perform 
repeated measurements under the same p requires the 
restoration of p before each measurement. 

In particular, we will consider a canonical ensemble 
of quantum systems at temperature 1/P (setting Boltz- 
mann's constant to 1) characterized by the density op- 
erator p — exp(— /3_ff)/(Tr exp(— Furthermore, we 
assume a Hamiltonian H = T + V being the sum of a 
kinetic energy term T = — (l/2m)A for a particle with 
mass m (A denoting the Laplacian and h = I) and an 
unknown local potential V{x,x') = v{x)S{x — x') which 
we want to reconstruct from data. In case of repeated 
measurements in a canonical ensemble one has to wait 
with the next measurement until thermal equilibrium is 
reached again. We will in the following focus on measure- 
ments of particle coordinates in a single particle system 
in a heath bath with temperature In that case the 
Xi represent particle positions corresponding to measure- 
ments of the observable Oi = x with x\xi> = Xi\xi>. 
The likelihood becomes 

p{xi\x,v) = '^Pa\(j)a{Xi)\^ =< \(f>{Xi)\'^ >, (3) 
a 

for H\(f>ci> = Ea\(t>a> and with <■ • •> denoting a thermal 
expectation for probabilities Pa — exp{—f3Ea)/Z with Z 
= exp(— One may also add a classical noise 
process p{xi\xi) on top of Xi in which case an additional 
integration is necessary to get the density p{xi\x,v) — 
Jdxip{xi\xi) p{xi\x,v) of the noisy output Xi. 

Even at zero temperature complete knowledge of the 
true likelihood would not be enough to determine a po- 
tential uniquely. The situation is still worse in prac- 
tice where only a finite number of probabilistic measure- 
ments is available and at finite temperatures, as the like- 
lihood becomes uniform in the infinite temperature limit. 
Hence, in addition to Eq.(|^) giving the likelihood model 
of quantum mechanics, it is essential to include a prior 
density over a space of possible potentials v. To be able 
to formulate a priori information explicitly, i.e., in terms 
of the functions values v{x) itself, we use a stochastic 
process. Technically convenient is for example a Gaus- 
sian process prior having the form 

j_ 

Po{v) = (det Ako") ' e-|<--o|Ko|.-.o>_ (4) 



with mean Vq, representing a reference potential or tem- 
plate for V, and a real symmetric, positive (semi-) definite 
covariance operator (AKq)^^, acting on potentials v and 
not on wave functions and measuring the distance be- 
tween V and vq. Hereby A is technically equivalent to 
a Tikhonov regularization parameter. Typical choices 
for Ko implementing smoothness priors are the negative 
Laplacian Kq = —A or a Radial Basis Function prior Kq 
= exp (— cT^gpA/2) More general, Gaussian process 
priors can be related to general approximate symmetries. 
Assume, for example, we expect the potential to com- 
mute approximately with a unitary symmetry operation 
S. Then, V ^ S'<VS ^ SV, defining an operator S act- 
ing on potentials V. In that case a natural prior would 
be po oc exp{-Es) with Es ^ <V - SV\V ~ SV > /2 
= <F|KoH/>/2 for Ko = (I-S)t(I-S) and I denot- 
ing the identity. Note that symmetric potentials are in 
the null space of such a Kq so also another prior has to 
be included if not already the combination with training 
data does determine the potential. Similarly, for a Lie 
group S{9) — exp (0s) an approximate infinitesimal sym- 
metry is implemented by Kq = s'^^s. In particular, a neg- 
ative Laplacian smoothness prior enforces approximate 
symmetry under infinitesimal translation. Alternatively, 
a more explicit prior implementing an approximate sym- 
metry can be obtained by choosing a symmetric reference 
potential Vs = SWsS and Es <V ~ Vs\V ~-Vs> /2. 

While a Gaussian process prior is only able to model 
a simple concave prior density, more general prior densi- 
ties can be arbitrarily well approximated by mixtures of 
Gaussian process priors 

M 

Po{v) ^^PkPo{v\k), (5) 
fc 

with mixture probabilities pk and Gaussian prior pro- 
cesses po('^lA) having means Vk and covariances (AK^)"^, 
where A plays the role of a "mixture temperature" . Simi- 
larly, one may introduce so called hyperparameters which 
parameterize the prior density [l^ . Those hyperparame- 
ters have to be included in the set of integration variables 
to obtain the predictive density and can also be treated 
in maximum a posteriori approximation. 

To find the potential with maximal posterior we have 
to maximize 

n 

p{v\D) OC po{v)Y[p{xt\Oi,v), (6) 

i 

where we assumed independent training data. This can 
be done by setting the functional derivatives of the log- 
posterior (technically often more convenient than the 
posterior) with respect to v{x) to zero, i.e., 

6,(^,^\npiv\D) = (7) 

where Sy(^x) stands for 6/ {Sv{x)). To calculate the func- 
tional derivative of the likelihoods in Eq.(Bp we require 
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5v(x)Ea as well as 5y(^x)<i)a{x'). Those can be found by 
taking the functional derivative of the eigenvalue equa- 
tion of the Hamiltonian yielding, for orthonormal eigen- 
functions, 

5v(x)Ea = <4>a\Sy(x)H \4>a>= |(/)a(a;)p, (8) 

Sv{x)(l}a{x') = ^ -—-^—(t)j{x')(l)*{x)(f>a{x), (9) 

using the fact that S^(^x)H{x' , x") — Sy(^^-^V{x' ,x") — 
5{x~x')5{x' — x"). Now it is straightforward to calculate 
the functional derivative of the likelihood 

5v{x)P{x^\x,v) = (10) 
< {5v(x)4'* {Xi)) H^i) > + < (j)* {Xi)S^,^^-)(j){x^) > 

-(3{< > - < |0(x,)r >< \4>\x)\'' >) . 

Finally, using 



Inpo = -AKo(?; - vq) 



(11) 



the functional derivative of the posterior can be calcu- 
lated. The formula ( p] ) is also valid for Gaussian mix- 
ture models of the form (^) provided we understand Kq 
= Z]fcPo(fc|«) Kfc and wq = K(7^ Po(fc|w) K^Wfe where 
Po{k\v) = Po{v\k)po{k)/pQ{v). 

It is straightforward to include also other kinds of 
data or a priori information. For example, a Gaussian 
smoothness prior as in Eq.(^ with zero reference poten- 
tial vo = and, say, v{x) = at the boundaries tends to 
lead to flat potentials when the regularization parameters 
A becomes large. For such cases it is useful to include be- 
sides smoothness also a priori information or data which 
are related more directly to the depth of the potential. 
One such possibility is to include information about the 
average energy U ^ < E > —J2aPaEa. The average 
energy can then be controled by introducing a Lagrange 
multiplier term ii{U — k), or, technically sometimes eas- 
ier, by a term representing noisy energy data 
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so that U K for /i — > oo. Using 5^(^x)U = < 6y(^x)E > 
~(3 <E Sy(^^)E> + P <E><Sy(^^)E>, we find for the 
functional derivative of Eu 



v{x)Eu = M 



(C/-k)(<|0(x)|2> (13) 

|0(x)|2>-t/<|0(x)|2>)). 

Collecting all terms we are now able to solve the sta- 
tionarity Eq.(|^) by iteration. Starting with an initial 
guess v^^\ choosing a step width ry and a positive defi- 
nite iteration matrix A we can iterate according to 



W + ,7A-i(Kf(i;o -!;('=)) (14) 
-^41np(x,|i,fW)-4£;[ 



Here we included an Ejj term which depends on v, like 
Kq for mixture models, and thus changes during itera- 
tion. Typically, r] and often also A are adapted during 
iteration. In the numerical examples we have studied A 
— XKq proved to be a good choice. 

The numerical difficulties of the nonparametric 
Bayesian approach arises from the fact that the quantum 
mechanical likelihood (^ is non-Gaussian and non-local 
in the potential v{x). Thus, even for Gaussian priors no 
v(a;)-integration can be carried out analytically. In con- 
trast, for example, Gaussian regression problems have a 
likelihood being Gaussian and local in the function of 
interest, and an analogous nonparametric Bayesian ap- 
proach with a Gaussian process prior requires only to deal 
with matrices with dimension not larger than the number 
of training data ||l^. The following examples will show, 
however, that a direct numerical solution of Eq.(|^) by 
discretization is feasible for one-dimensional problems. 
Higher dimensional problems, on the other hand, require 
further approximations. For example, treating inverse 
many-body problems on the basis of a Hartree-Fock ap- 
proximation is ongoing work. 



(a) likelihoods 



(b) potentials 





FIG. 1. Reconstruction of an approximately periodic po- 
tential from coordinate measurements, (a) shows empirical 
density (bars), true likelihood (thin line), and reconstructed 
likelihood (thick line), (b) shows the true potential (thin line), 
reference potential vo = sin(7ra::/3) (dashed line), and recon- 
structed potential v (thick line). (200 data points, mesh with 
30 points, m = 0.25, f3 = A, Kq = -A, A = 0.2, /i = 0, 
(7(wtruc) = -0.354 resuhing in U{v) = -0.552.) 



(a) potentials 



(b) potentials 





FIG. 2. Reconstruction with energy penalty term Eu with 
K = (/(titruo) = —0.354. Shown are true potential (thin 
lines), reconstructed potential (thick lines), reference poten- 
tial (dashed), (a) Ko ~ —A, vo{x) — sin(7r2;/3), fj, = 1000, 
U{v) = -0.353. (b) Ko = -A - -yAr, vq = 0, T = 6, -y = 
0.12, A = 0.05, = 10, U{v) = -0.351. All other parameters 
as for Fig.|l| 



3 



(a) likelihoods (b) potentials 




FIG. 3. Reconstruction of symmetric potential with mix- 
ture of Gaussian priors, (a) shows empirical density (bars), 
true likelihood (thin line) , and reconstructed likelihood (thick 
line), (b) shows the true potential (thick line), the two refer- 
ence potential vi and V2 (dashed lines), and the reconstructed 
potential v (thick line). (20 data points, mesh with 31 points, 
for v{x) = v{—x) and boundary condition w(±15) = 0, with m 
= 0.1, P = l,Ko = -A, A = 0.5, ^ = 10 «: = -9.66 = (7(«truo) 
resulting in U{v) = —10.1. ) Because both reference poten- 
tials are partially supported by the data the approximated v 
is essentially a smoothed mixture between vi and «2 in regions 
where no data are available with mixture coefficients for prior 
components po{l\v) — 0.57, po{2\v) = 0.43. 

In the following numerical examples we discuss 
the reconstruction of an approximately periodic, one- 
dimensional potential. As a first example, assume we 
expect V to be similar to a periodic reference potential 
vo, which we choose as reference potential. To enforce 
the deviation from vo to be smooth we take as prior on 
V a negative Laplacian covariance, i.e., Kq = —A. Fig.|l| 
shows representative numerical results for a grid with 
30 points and 200 data sampled from a true likelihood 
p{x\x , Vfmc) corresponding to a true potential z;true. The 
reconstructed potential v have been found by iterating 
without energy penalty term Eij according to Eq.(p^) 
with A = Kq. Choosing zero boundary conditions for v 
the matrix Kg is invertible. Consistent with the bound- 
ary conditions on v we took periodic boundary conditions 
for the eigenfunctions Note that the data have been 
sufhcient to identify clearly the deviation from the peri- 
odic reference potential. Fig.^ shows the same exam- 
ple with an energy penalty term Eu with k = C/(ftruo)- 
While the reconstructed likelihood (not shown) is not 
much altered, the true potential is now better approxi- 
mated in regions where it is small. 

Fig.^ shows the implementation of approximate pe- 
riodicity by an operator At measuring the difference 
between the potential v{x) and the potential trans- 
lated by T and and which is defined by <w|A7i|i;> — 
Jdx\v{x)~v{x+T)\'^ for periodic boundary conditions on 
V. To find smooth solutions we added a negative Lapla- 
cian term with zero reference potential, i.e., we used Kq 
= — A — 7 At- To have an invertible matrix for periodic 
boundary conditions on v we iterated this time with A = 
Kq + O.ll. The implementation of approximate period- 
icity by At instead of an periodic vq is more general in as 
far it allows arbitrary functions with period T. As, how- 
ever, the reference function of the Laplacian term does 
not fit the true potential very well, the reconstruction is 



poorer in regions where the potential is large and thus no 
data are available. In these regions a priori information 
is of special importance. Finally, Fig.^ shows the imple- 
mentation of a mixture of Gaussian prior processes. 

In conclusion, we have applied a nonparametric 
Bayesian approach to inverse quantum statistics and 
shown its numerical feasibility for one-dimensional ex- 
amples. 
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